COMMENT ⊗   VALID 00008 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	% This program by D. E. Knuth is not copyrighted and can be used freely.
C00004 00003	@* Introduction.
C00011 00004	@* Inputting and outputting the data.
C00023 00005	@* Error diffusion.
C00035 00006	@* Computing the diffusion tables.
C00041 00007	@* The main program.
C00043 00008	@* Index.
C00054 ENDMK
C⊗;
% This program by D. E. Knuth is not copyrighted and can be used freely.

% Here is TeX material that gets inserted after \input webmac
\def\title{PICFONT}
\font\logo=logo10
\def\MF{{\logo META}\-{\logo FONT}}

%\advance\topskip by \baselineskip	% doublespacing
%\advance\smallskipamount by \baselineskip
%\advance\baselineskip by \baselineskip

\def\con{\par\vfill\eject % finish the section names
  \rightskip 0pt \hyphenpenalty 50 \tolerance 200
  \setpage
  \output{\normaloutput\page\lheader\rheader}
  \titletrue % prepare to output the table of contents
  \pageno=\contentspagenumber \def\rhead{TABLE OF CONTENTS}
  \message{Table of contents:}
  \topofcontents
  \line{{\bf Sample}\hfil Section}
  \def\Z##1##2##3{\line{\ignorespaces##1
    \leaders\hbox to .5em{.\hfil}\hfil\hbox to2em{\hss##2}}}
  \readcontents\relax % read the contents info
  \botofcontents \end} % print the contents page(s) and terminate
@* Introduction.
This program prepares a \MF\ program for a special-purpose font
that will approximate a given picture.
The input is assumed to be a binary file that contains one byte of
density information per pixel.
The output will be a sequence of lines like
$$\hbox{\tt row(10); cols(3,15,16,17);}$$
this means bits 3, 15, 16, and 17 of the character for row 10 should be black.

@ Here's an outline of the entire Pascal program:

@p program picfont(@!bytes_in,@!output);
type @<Types in the outer block@>@/
var@?@<Global variables@>@/
@<Basic procedures@>
begin @<The main program@>;
end.

@ The picture in the input data is assumed to contain |mm| rows and |nn| columns.

@d mm=512 {this many rows}
@d nn=440 {this many columns}

@ It's convenient to declare a macro for incrementation.

@d incr(#) == #←#+1
@* Inputting and outputting the data.
The input appears in a file of 8-bit bytes, with \.{00} representing black
and \.{FF} representing white. There are $mm\times nn$ bytes; they appear in
order from top to bottom and left to right just as we normally read a page
of text.

@<Types...@>=
@!eight_bits=0..255; {unsigned one-byte quantity}
@!byte_file=packed file of eight_bits; {files that contain binary data}

@ @<Glob...@>=
@!bytes_in:byte_file;

@ Different Pascal systems have different ways of dealing with
binary files. Here is one common way.
@↑system dependencies@>

@<Open the input file@>=
reset(bytes_in,'','/B:8')

@ We shall use the following model for estimating the effect of a
given bit pattern: If a pixel is black, the darkness is 1.0; if it
is white but at least one of its four neighbors is black, the darkness
is |zeta|; if it is white and has four white neighbors, the darkness
is zero.

@d white=0 {code for a white pixel with all white neighbors}
@d gray=1 {code for a white pixel with 1, 2, 3, or 4 black neighbors}
@d black=2 {code for a black pixel}
@d zeta==0.2 {assumed darkness of white pixel with a black neighbor}

@ There isn't room to store all the input bytes in memory at once, but
it suffices to keep buffers for about a dozen rows near the current area
being computed.

@<Glob...@>=
@!ii:integer; {the buffer holds rows |8ii-7| through |8ii+4|}
@!buffer:array[-2..9,0..nn+1] of real; {densities in twelve current rows}
@!darkness:array[-3..9,0..nn+1] of white..black; {darknesses in buffer rows}
@!new_row:array[0..nn+1] of real; {densities in row being input}

@ The |get_in| procedure computes the densities in a specified row
and puts them in |new_row|. This procedure is called successively for
|i=1|, 2,~\dots\thinspace.

@<Basic procedures@>=
procedure get_in(@!i:integer);
var @!j:integer;
@!t:eight_bits; {byte of input}
begin new_row[0]←0.0;
if i>mm then for j←1 to nn do new_row[j]←0.0
else for j←1 to nn do
	begin read(bytes_in,t); new_row[j]←(255.5-t)/256.0;
	end;
new_row[nn+1]←0.0;
end;

@ Here is a procedure that ``rolls'' the buffer down eight lines:

@<Basic procedures@>=
procedure roll;
var @!j:0..nn+1;
@!i:2..9;
@!k:integer;
begin for i←6 to 9 do for j←0 to nn+1 do
	begin buffer[i-8,j]←buffer[i,j]; darkness[i-8,j]←darkness[i,j];
	end;
for j←0 to nn+1 do darkness[-3,j]←darkness[5,j];
incr(ii);
for i←2 to 9 do
	begin get_in(8*ii+i-5);
	for j←0 to nn+1 do
		begin buffer[i,j]←new_row[j]; darkness[i,j]←white;
		end;
	end;
end;

@ It's tedious but not difficult to get everything started.
We put zeros above the top lines in the picture.

@<Initialize the buffers@>=
ii←0;
for i←6 to 9 do
	begin get_in(i-5);
	for j←0 to nn+1 do
		begin buffer[i,j]←new_row[j]; darkness[i,j]←white;
		end;
	end;
for i←-2 to 5 do for j←0 to nn+1 do
	begin buffer[i,j]←0.0; darkness[i,j]←white;
	end;
for j←0 to nn+1 do darkness[-3,j]←white

@ It's easy to output the current darkness values. Here we output
eight consecutive rows.

@<Output the pixel values for the top eight rows of the buffer@>=
for i←-2 to 5 do
	begin write('row(',8*ii-5+i:1,'); cols('); cols_out←0;
	for j←1 to nn do if darkness[i,j]=black then
		begin if cols_out<15 then
			begin if cols_out>0 then write(',');
			incr(cols_out);
			end
		else	begin write_ln(','); write('               '); cols_out←1;
			end;
		write(j:1);
		end;
	write_ln(');')
	end

@ @<Glob...@>=
@!cols_out:0..15; {the number of columns output so far on this line}
@* Error diffusion.
The pixels are divided into 64 classes, numbered from 0 to~63. 
We convert the pixel values to 
darknesses by assigning them for class~0 first, then class~1,
etc. The error incurred at each step is distributed to the neighbors whose
class numbers are higher. This is done by means of precomputed tables
|class_row|, |class_col|, |start|, |del_i|, |del_j|, and |alpha| whose
function is easy to deduce from the following code:

@<Choose pixel values and diffuse the errors in the buffer@>=
for k←0 to 63 do
	begin i←class_row[k]; j←class_col[k];
	while j≤nn do
		begin @<Decide the color of pixel |[i,j]| and the resulting |err|@>;
		for l←start[k] to start[k+1]-1 do
			begin u←i+del_i[l]; v←j+del_j[l];
			buffer[u,v]←buffer[u,v]+err*alpha[l];
			end;
		j←j+8;
		end;
	end

@ @<Glob...@>=
@!class_row:array[0..63] of -2..8;
	{buffer row containing pixels of a given class}
@!class_col:array[0..63] of 1..8;
	{first column containing pixels of a given class}
@!class_number:array[-2..9,0..9] of 0..63; {number of a given position}
@!err:real; {error introduced at the current position}
@!err_black:real; {error introduced at the current position if black chosen}
@!black_diff:real; {difference between |err| and |err_black| for gray pixel}
@!l:0..256; {index into diffusion tables}
@!start:array[0..64] of 0..256;
	{first entry of diffusion table for a given class}
@!del_i,@!del_j:array[0..256] of -1..1; {neighboring location for diffusion}
@!alpha:array[0..256] of real; {constant of proportionality for diffusion}

@ Here we choose white or black, whichever minimizes the magnitude of the error.
The |gray| values of this pixel and its neighbors make this calculation slightly
tricky, as we must subtract |zeta| when a gray pixel is created and add |zeta|
when it is destroyed.

@<Decide the color of pixel |[i,j]| and the resulting |err|@>=
if darkness[i,j]=gray then
	begin err←buffer[i,j]-zeta; err_black←err-black_diff;
	end
else	begin err←buffer[i,j]; err_black←err-1.0;
	end;
if darkness[i-1,j]=white then err_black←err_black-zeta;
if darkness[i,j-1]=white then err_black←err_black-zeta;
if darkness[i,j+1]=white then err_black←err_black-zeta;
if darkness[i+1,j]=white then err_black←err_black-zeta;
if err_black+err>0 then
	begin err←err_black; darkness[i,j]←black;
	if darkness[i-1,j]=white then darkness[i-1,j]←gray;
	if darkness[i,j-1]=white then darkness[i,j-1]←gray;
	if darkness[i,j+1]=white then darkness[i,j+1]←gray;
	if darkness[i+1,j]=white then darkness[i+1,j]←gray;
	end

@ @<Initialize the diffusion tables@>=
black_diff←1.0-2.0*zeta;
@* Computing the diffusion tables.
The tables for dot diffusion could be specified by a large number
of boring assignment statements, but it is more fun to compute them by a method
that shows some of the mysterious underlying structure.

@<Initialize the diffusion tables@>=
@<Initialize the class number matrix@>;
@<Compile ``instructions'' for the diffusion operations@>

@ The order of classes
used here is the order in which pixels might be blackened in a font
for halftones based on dots in a 45$↑\circ$ grid.

@<Basic procedures@>=
procedure store(@!i,@!j:integer); {establish new |class_row|, |class_col|}
begin if i<1 then i←i+8@+else if i>8 then i←i-8;
if j<1 then j←j+8@+else if j>8 then j←j-8;
class_number[i,j]←k; class_row[k]←i; class_col[k]←j; incr(k);
end;
@#
procedure store_eight(@!i,@!j:integer); {rotate and shift for eight classes}
begin store(i,j); store(i-4,j+4); store(5-j,i); store(1-j,i-4);@/
store(4+j,1-i); store(j,5-i); store(5-i,5-j); store(1-i,1-j);
end;

@ @<Initialize the class number matrix@>=
k←0; store_eight(7,2); store_eight(8,3); store_eight(8,2); store_eight(8,1);@/
store_eight(1,4); store_eight(1,3); store_eight(1,2); store_eight(2,3);@/
for i←1 to 8 do
	begin class_number[i,0]←class_number[i,8];
	class_number[i,9]←class_number[i,1];
	end;
for j←0 to 9 do
	begin class_number[-2,j]←class_number[6,j];
	class_number[-1,j]←class_number[7,j];
	class_number[0,j]←class_number[8,j];
	class_number[9,j]←class_number[1,j];
	end

@ The tricky part of this process is the fact that some values near the
bottom of the buffer aren't ready for processing until errors have been
diffused from the next bufferload. In such cases we go up eight rows
to process a value that has been held over.

@<Glob...@>=
@!hold:array[0..9,0..9] of boolean; {is this value too close to the bottom
	of the buffer to allow immediate processing?}

@ The ``compilation'' in this step simulates going through the diffusion
process the slow way, and records the actions it does (so that they
can all be done a high speed later).

@<Compile...@>=
for j←0 to 9 do hold[9,j]←true;
for i←0 to 8 do for j←0 to 9 do hold[i,j]←false;
l←0; k←0;
repeat i←class_row[k]; j←class_col[k]; w←0; start[k]←l;
for u←i-1 to i+1 do for v←j-1 to j+1 do
 if class_number[u,v]>k then
	begin del_i[l]←u-i; del_j[l]←v-j; incr(l);
	if u=i then w←w+2 {neighbors in the same row get weight 2}
	else if v=j then w←w+2 {neighbors in the same column get weight 2}
	else incr(w); {diagonal neighbors get weight 1}
	end
 else	if hold[u,v] then hold[i,j]←true;
if hold[i,j] then class_row[k]←i-8;
@<Compute the |alpha| values for class |k|, given the total weight |w|@>;
incr(k);
until k=64;
start[64]←l

@ @<Compute the |alpha| values for class |k|, given the total weight |w|@>=
for ll←start[k] to l-1 do
	begin if del_i[ll]=0 then alpha[ll]←2.0/w
	else if del_j[ll]=0 then alpha[ll]←2.0/w
	else alpha[ll]←1.0/w;
	end

@ @<Glob...@>=
@!ll:0..256; {loop index}
@!u,@!v:integer; {neighbors of |i| and |j|}
@!w:integer; {the weighted number of high-class neighbors}
@!i,@!j:integer; {the current pixel position being considered}
@!k:0..64; {the current class being considered}
@* The main program.
Finally we're ready to get it all together.

@<The main program@>=
@<Initialize the diffusion tables@>;
@<Open the input file@>;
@<Initialize the buffers@>;
repeat @<Choose pixel values and diffuse the errors in the buffer@>;
if ii>0 then @<Output the pix...@>;
roll;
until 8*ii>mm
@* Index.
Here are the quantities declared and/or used in the program.
(The uses of single-letter variables aren't indexed.)